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ABSTRACT 

Black holes that accrete far below the Eddington limit are believed to accrete through a geometrically 
thick, optically thin, rotationally supported plasma that we will refer to as a radiatively inefficient 
accretion flow (RIAE). RIAEs are typically collisionless in the sense that the Coulomb mean free path 
is large compared to GM/c^, and relativistically hot near the event horizon. In this paper we develop 
a phenomenological model for the plasma in RIAEs, motivated by the application to sources such as 
Sgr A* and M87. The model is derived using Israel-Stewart theory, which considers deviations up to 
second order from thermal equilibrium, but modified for a magnetized plasma. This leads to thermal 
conduction along magnetic field lines and a difference in pressure, parallel and perpendicular to the 
field lines (which is equivalent to anisotrotropic viscosity). In the non-relativistic limit, our model 
reduces to the widely used Braginskii theory of magnetized, weakly collisional plasmas. We compare 
our model to the existing literature on dissipative relativistic fluids, describe the linear theory of the 
plasma, and elucidate the physical meaning of the free parameters in the model. We also describe 
limits of the model when the conduction is saturated and when the viscosity implies a large pressure 
anisotropy. In future work, the formalism developed in this paper will be used in numerical models 
of RIAEs to assess the importance of non-ideal processes for the dynamics and radiative properties of 
slowly accreting black holes. 


1. INTRODUCTION AND ASTROPHYSICAL 
CONTEXT 

Most massive galaxies have black holes at their cen¬ 
ters, and most of these black holes are accreting far 
below the Eddington rate MEdd Low lu¬ 

minosity black holes are believed to accrete through a 
geometrically thick, optically thin disk. Phenomeno¬ 
logical (r adiatively inefficient a ccretion flows or RI- 
AEs, see lYuan fc NaravanI 1201^ and numerical (gen¬ 
eral relativ istic magnetohydrodynamics, or GRMHD, 
see e.g. iKoide et alJ II999I : iDe Villiers et al.l 120031 : 
iMcKinnev fc Cammid I2004H models suggest that the 
density and temperature of the accreting plasma are 
such that the collisional (Coulomb) mean free paths (ion- 
ion, ion-electron, and electron-electron) are many orders 
of magnitude larger than GM /c^ when M <g; MEdd 
(|Mahadevan fc Quataert! [1 9 9711 . The accreting plasma 
is thus collisionless. 

In the nearby universe, the roster of low luminosity 

manic@illinois.edu 
gammie@illinois.edu 
tvtoucart (Slbi.gov 
eliot @ berkeiey. edu 


black holes includes M87 (accreting 4-5 orders of magni¬ 
tude below MEdd), and Sgr A* (accreting about 8 orders 
of magnitude below MEdd)- These two sources are the 
largest known black holes in terms of angular size on 
the sky. As a result, they are the two main targets for 
high resolut ion imaging experiment s, including Gravity 
on the VLT (|Eisenha uer et a l]l2008f) an d the Event Hori¬ 
zon Telescope (EHT: iDoeleman et al.l 120091) : the latter 
will use submillimeter very long baseline interferometry 
(VLBI) to resolve the accretion flow and jet on angular 
scales comparable to the event horizon. The EHT data 
may be used to test General Relativity by measuring the 
angular size of the photon orbit, if astrophys ical uncer¬ 
tainties can be controlled (iPsaltis et al.ll20Lll . 

The fact that the collisional mean free path is much 
larger than GM/c^ in RIAEs implies that non-ideal pro¬ 
cesses such as conduction and viscosity are likely to be 
important. Furthermore, the mean free path is also much 
larger than the Larmor radii of all the species in the 
plasma and the gyration time scale is much shorter than 
the dynamical time scale. This leads to the above dis¬ 
sipative processes being anisotropic with respect to the 
local magnetic field. Heat flows only along the field lines 
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and an anisotropic viscosity is generated by a shear flow 
projected along the field lines. 

Consider conduction, using non-relativistic estimates: 
the volume heating rate is du/dt = — V-q, where q is the 
conductive heat flux. The timescale for changing the in¬ 
ternal energy is Tcond = u/{du/dt). Estimate V-q ^ g/r, 
assume the heat flux is approximately saturated (given 
the large mean free paths) so that q ^ uv, where u is 
the internal energy density and v is rms particle speed. 
Then Tcond ^ t/c for electrons (electrons are relativistic 
close to the horizon, so z; ~ c) and Tcond ^ r/cs for pro¬ 
tons, which in a near-virial RIAF is the dynamical time. 
Conduction can thus potentially play an important role 
in co ntrolling the thermal state of the accreting plasma 
(e.g., [Johnson fc Quataertll2007f) . 

What about viscosity or equivalently (as we show be¬ 
low) a difference in pressures along and perpendicular to 
the local magnetic field? Pressure anisotropy can be gen¬ 
erated in a number of ways, for example, by anisotropic 
compression or expansion of the plasma, through a linear 
shear, etc. If the magnetic field strength varies, adiabatic 
invariance of the magnetic moment associated with the 
orbit of charged particles about a magnetic field with 
strength B implies that T±/B is invariant for a plasma 
with kT < mc^ {T± = temperature perpendicular to the 
local magnetic field). Here, the temperature T and the 
mass m corresponds to a specific species. Thus, if the 
plasma is compressed in the plane perpendicular to the 
mean field so that the density increases by a factor R, 
the perpendicular temperature also increases by a fac¬ 
tor R, generating a significant pressure anisotropy, i.e., 
viscosity. Order unity fluctuations in density and mag¬ 
netic field strength are common in numerical models of 
accretion flows, so order unity pressure anisotropy is ex¬ 
pected in the absence of collective effects. This implies 
that viscous stresses may be dynamically important and 
contribute significan tly to angular momentum transport 
and plasma heating (jSharma et al.l[2?i0^ I2007D . 

Despite their potential importance, conduction and 
viscosity are, however, absent from all global relativistic 
numerical models of RIAFs to date. This is one of the 
significant systematic uncertainties in developing models 
for the emission from systems such as Sgr A* and M87. 

In this paper, we develop a formalism for modeling rel¬ 
ativistic anisotropic conduction and viscosity, motivated 
by the application to RIAFs. Although the plasmas of 
interest are macroscopically collisionless, we focus in this 
paper on the more modest task of developing a theory for 
collisional magnetized plasmas in which the mean free 
path is large compared to the Farmer radius of particles, 
but small compared to the system size. The former hi¬ 
erarchy implies that heat and momentum transport are 
predominantly along the local magnetic field direction, 
while the latter constraint implies that one can derive 
the relevant equations using an expansion about thermal 
equilibrium. Our assu med hierarchy of length-scales is 
similar to that used in lBraginskiillI96^ s theory of non- 
relativistic magnetized plasmas, which has been widely 
applied to ui iderstand the p hysics of dilute astrophysical 
plasmas Isee lKulsrudll2005D . We will refer to our formal¬ 
ism as an extended MHD (EMHD) model. 

Although the applications that motivate this work are 
to collisionless systems, we focus on the collisional regime 


for the following reasons: (I) the th eory of dissipative rel¬ 
ativis tic fluids is quite subtle (e.g. lAndersson fc Comer! 
[2001 . so it seems prudent to not jump directly to the 
yet more challenging long mean free path regime; (2) 
wave-particle interactions and velocity space instabilities 
limit the mean free path of charged particles to be much 
less than the colli sional mean free path under the con¬ 
ditions of interest (iSharma et al.ll^006HKunz et al.ll20l4 
iRiauelme et al.ir20I5[l implying that the ‘collisional’ the¬ 

ory may be more appropriate than one might have first 
anticipated. The formalism we develop allows the vis¬ 
cosity and conductivity to depend arbitrarily on local 
plasma conditions so that these wave-particle limits on 
the mean free path can be incorporated as sub-grid mod¬ 
els. 

Throughout the paper we formulate the equations 
in terms of a single fluid model. In reality, the low- 
collisionality plasmas of interest are believed to develop 
a two-temperature structure because the timescale for 
Coulomb collisions to equilibrate the electron and proton 
temperatures is long compared to the dynamical time in 
the accreting plasma. A formulation for dealing with 
electron dynamics and its n umerical implementa tion has 
been recently introduced bv iRessler et ah! (j2015[ l. where 
a reduced form of our conduction model has been used 
and appropriately modified for electrons. 

The remainder of this paper is organized as follows. 
In §2 we write down basic equations and describe the 
equivalence of viscosity and anisotropic pressure. In §3 
we describe the desired asymptotic behavior of any ex¬ 
tended MHD closure model. In §4 we derive the evolu¬ 
tion equations for the heat flux and pressure anisotropy. 
§5 describes the connection between our model and non- 
relativistic dissipative theory. §6 motivates a scheme for 
fixing the model parameters (the transport coefficients : 
viscosity and conductivity) in terms of a relaxation time. 
§7 gives the linear theory and the stability thresholds of 
our model, and §8 a brief discussion of nonlinear (shock) 
solutions. Finally in §9 we offer a guide to the model, a 
summary of the formalism, and the relationship to earlier 
works. In the Appendix we show how the characteristic 
pressure anisotropy derived on thermodynamic grounds 
using the Israel-Stewart theory can also be interpreted 
as arising from conservation of relativistic adiabatic in¬ 
variants in a magnetized plasma. 

2. PHYSICAL CONTEXT 

We begin by defining notation and frames. 

We work in a spacetime described by the metric 
whose determinant we denote with g. Consider a plasma 
consisting of particles with distribution function fs = 
dN/d^xd^p, and rest mass mg, where d^p = dpidp 2 dp 3 , 
Pi are the spatial covariant components of the particle 
four-momentum and s indicates the species of the 
particles (electrons, ions, etc.). The distribution function 
is invariant. Each species has a number current 

We assume that the plasma consists of electrons and ions, 
and is quasi-neutral everywhere. Thus both these species 
have the same approximate number current = Nj/ « 
Nj/. We define the rest frame as that in which the num- 
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her current has no spatial components, 
have, 


u 


M = 




n 


Therefore we 


( 2 ) 


where n is the number density of ions, and is equal to the 
number density of electrons. This definition of implies 
that we are using the so-called Eckart frame, in which 
mass diffusion is absent. An alternative, the Landau 
frame, assumes that energy diffusion is absent. The total 
rest mass density p, often denoted as po in the relativity 
literature is p = —TismsNj^u^ = —{rmN^ 

Since « iV^ = we have p = —{rm -\- me)N^^u^ = 
(rrii -I- me)n. 

The matter stress-energy tensor is 


^ matter — 



d^p 


P^P^fs- 


( 3 ) 


In the fluid frame, the spatial components of the ideal 
fluid pressure tensor are 


/ P 0 0 \ 

pb = 0 P 0 . (12) 

V 0 0 P / 

In what follows we are interested in modeling a magne¬ 
tized plasma that departs from ideality in that it has 
a conductive heat flux and a viscous stress. The to¬ 
tal stress tensor with the fluid assumed perfect and 
the electromagnetic terms included under the ideal mag¬ 
netohydrodynamics (MHD) approximation (conductivity 
a = oo) is 

= (p + w + ^b^)u^u'' + (P + y- 6^6". (13) 
where 6^ = 6^6^ and 


Each of these definitions is invariant because 
d^p/{y/—gp*) is invariant. 

On taking moments of the Boltzmann equation, one 
can show that the quantities and Tj^atter satisfy 

= V^(nu'^) = 0 (4) 

and 

^^TZtter = ( 5 ) 

where is the covariant derivative, is the electro¬ 
magnetic field tensor and is the four-current. The 
divergence of the electromagnetic stress tensor is 

( 6 ) 

On adding dS]) and ([6]), we get that the divergence of the 
total stress tensor is zero, as required by the Bianchi 
identities 


where e = Levi-Civita tensor, which is antisymmetric on 
all pairs of indices. Evidently = 0, and b^^ reduces 
to the magnetic field in the fluid frame (with a factor of 
•\/4^ absorbed into the definition). The magnetic field 
evolution is given by 

- 6'^m") = 0 (15) 

which combines the induction equation (three space com¬ 
ponents) with the no-monopoles condition (time compo¬ 
nent). 

When including non-ideal effects in the stress-energy 
tensor, the heat flux in the system is 

= -KupT^^. (16) 


= = ^ ( 7 ) 

For a perfect unmagnetized fluid the stress-energy ten¬ 
sor is 

TZ.,,,r = {P + y)u^u'' + Ph>^'^, (8) 

where 

= gt^'' + u>^u'' (9) 


Combined with ©, this provides a kinetic theory defi¬ 
nition for the heat flux. Evidently u^q^ = 0, so in the 
fluid frame, = 0. From m one can show that the 
heat flux makes a contribution to the stress tensor of the 
form 

= (17) 

which in the fluid frame has the form 


is the projection tensor (projects into the space normal to 
the fluid four-velocity), is the metric, u the internal 

energy per unit proper volume, and P the gas pressure. 

We define the fluid frame as an orthonormal tetrad 
with time component and three additional 

spacelike basis vectors and e^^y In the fluid 

frame the stress-energy tensor is 


/p-t-wO 0 0\ 

rp{,a)(h) _ I 0 poo] 

matter Q 0 P 0 ’ 

Vo OOP/ 


( 10 ) 


which, with ©, provides a kinetic theory definition for 
the pressure and internal energy. 

The space-space part of the stress-energy tensor is the 
pressure, or stress, tensor 




/ 0 q^ qy qZ 

rp{a){b) _ I 0 0 0 ] 

^cond - 0 0 0 

V 0 0 0 / 


(18) 


The viscous stress tensor models momentum fluxes 
set up by departures from equilibrium due to a shear flow. 
It is given by 


(19) 


It is perhaps not as widely appreciated as it should 
be, that the viscous stress can be recast as a pressure 
anisotropy. In the fluid frame the viscous stress tensor 
is a symmetric matrix, so there is always a basis (ob¬ 
tained by rotation) where it can be written in diagonal 
form: 


0 0 
Py 0 
OOP, 


0 


( 11 ) 


plj — 


( 20 ) 
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where there is a separate pressure for each direction. Be¬ 
low we consider a magnetized plasma where the stress 
tensor is symmetric under rotations around the mag- 


J 


1 , 


netic held. If the held is in the 2 ; direction, this implies 


Prr = Py 


With conduction and viscosity of the plasma included, 
the total stress-energy tensor is now 


1 , 


T^^Lter+EM =T^^’' = {p + u+ + (P + + q’^ 


( 21 ) 


This stress-energy tensor is quite general, but one needs 
an appropriate model for u, P, q^, and Note that 
the electromagnetic terms are still written down in the 
inhnite conductivity limit. 


3. MODEL DESIDERATA 


What are the desirable properties of a closure model 
for the heat hux q^ and the viscous stress 11^'^? 

(1) The model should be causal. For non-relativistic 
shear viscosity and thermal conduction, the energy and 
momentum huxes are proportional to gradients of the 
temperature and velocity, and so respond instanta¬ 
neously to changes in the huid. Classical, non-relativistic 
models are parabolic and have characteristics that propa¬ 
gate at inhnite speed. The classical theories ca n be made 
causal in a model pioneered by Maxwell and ICattane^ 
(|I948f) in which the energy and momentum huxes relax 
to their classical values on a characteristic timescale r. 

(2) The model should be stable. The relativistic ther¬ 
mal conduction model of lEckartI (119401) (see also MTW) 
S6tS 

qf^ = -pxh^^^{^,e + ea,) ( 22 ) 


where x is the thermal diffusivity, 0 = kTImc? 

is the normalized temperature and ai, is the four- 
acceleration. The term proportional to the four- 
acceleration drives the temperature toward a con¬ 
stant redshifted temperature rather than a con¬ 
stant local temperature - a desirable effect - but it 
makes the theory unstable (iGarcia-Perciante et al.ll2009t 

iLopez-Monsalvo fc AnderssonI 120111) . Long wavelength 
modes {k 0) are unstable with growth rate 


{pc^ -I- M -I- P)c? 

xP 


(23) 


Notice that as y —0, w —oo. If this theory were 
correct, the water in our bodies would spontan eously 
explode in 10“^"^ sec (|Hiscock fc LindblomI 1198^ . Ev¬ 
idently the stability of relativistic conduction theories 
is nontrivial. A relativistic extension of the Maxwell- 
Cattaneo procedure not only makes the theory hyper¬ 
bolic but also conditionally eliminates the Eckart insta¬ 
bility. 

(3) Entropy should increase, i.e. the model should obey 
the second law of thermodynamics. The entropy con¬ 
straint is expressed by defining an entropy four-current 
and requiring that > 0. This constraint was used 
to derive the Eckart model wherein the entropy current 
is expanded around equilibrium to first order in the heat 
flux. The first order model suffers from the instability 
de scribed above. Expandi ng up to second order, as done 
by [Israel fc Steward (119790 leads to conditionally hyper¬ 
bolic, stable and causal equations. We will use this in the 


next section to derive evolution equations for our model 
of anisotropic thermal conduction and viscosity. 

(4) We are interested in plasmas with ion and elec¬ 
tron Larmor radii tiny compared to the characteris¬ 
tic scale GM/c^, and ion and electron gyroperiod tiny 
compared to the dynamical timescale. Therefore, we 
shall assume that the distribution functions of both the 
ions and electrons are independent of the gyrophase, i.e. 
fs = fs{pi\,p±), where s indicates the species. Now in a 

tetrad frame with and = b^^, we eval¬ 

uate ([3|) to find only the following non-zero components 


rpiGG) _ 

matter 


/ y(0)(0) 

0 

0 

T°ll 


0 0 r°ii \ 

r-L-L 0 0 

0 t-l-l 0 

0 0 Tim / 


(24) 


The terms which are zero are identically so, because they 
appear as sm{6)dd or cos{6)dd, where 9 is the 
gyrophase. We see that there is a heat flux T^’H = q only 
along the magnetic field line. Therefore our model for 
the heat flux can be written as 


qi^=qbf^ (25) 

where b^ = b^ j y^Wb^. q will be the fundamental vari¬ 
able describing the heat fluxQ 
We now write down the pressure tensor with T±± = P± 
and Tim = T|| 



0 

0 ^ 

, / T + APq. 

0 

0 

pwU) — 0 

p± 

0 

= 0 

P-f APj_ 

0 

V 0 

0 

P\\J 

' [ 0 

0 

P + AP|| 


(26) 

where P is the ideal fluid pressure and AP± and APy 
are deviations from it in the T and || directions respec¬ 
tively. The variables APj^ and AP|| can in principle 
vary independently and give rise to both a bulk viscosity 
(trace) and a shear viscosity (trace-free part). We sim¬ 
plify the model further by assuming that bulk viscosity 
is zero and thus imposing that the deviation of the pres¬ 
sure tensor from ideality be trace-free. Doing so gives 
AP|| = —2AP^. Now redefining AP± = AP/3, we have 

^ A more accurate description of a collisionless plasma requires 
us to differentiate between a heat flow due to parallel temperatures 
gradients gji* and a heat flow due to perpendicular temperatures 

gradients (j*(, both of which flow along the field lines (/(( = g|| 6'' 

and q'^ = q^be. The net heat flow is then q>^ = qb^ = (q±+q\\)be. 
However, even in this case, the heat flux appears in the stress tensor 
only as the sum q = q± + qp 
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for the pressure tensor 

/ P+^AP 0 0 \ 

= 0 P+iAP 0 . (27) 

\ 0 0 P- |AP J 

We see that P± = P + AP/3 and P|| = P — 2AP/3. 
Therefore AP = P± — P||, which is the usual definition 
of pressure anisotropy. The shear stress in an arbitrary 
frame is then 


= - AP 



(28) 


AP is the fundamental variable describing the viscous 
stress. The above expression satisfies II)) = 0 and is thus 
trace-free. 

(5) If possible, the model should asymptote to a rig¬ 
orous model in the collisional limit. The relation to ear¬ 
lier theories will be discu ssed in detail late r , but in brief 
our model is equivalent to [Israel fc Stewa^ (II979I1 theory 
projected along the magnetic field lines. Israel-Stewart 
theory has 9 fields that describe nonideal effects: 5 for 
the shear viscosity, I for bulk viscosity, and 3 for the 
conductivity. Projecting along the magnetic field lines 
reduces the viscous shear stress degrees of freedom from 
5 to I, and the heat flux degrees of freedom from 3 to I, 
while we ignore bulk viscosity. 


4. EVOLUTION OF Q AND AP 

Following Israel and Stewart, it is possible to derive 
evolution equations for q and AP from the second law of 
thermodynamics, expressed here by the requirement that 
the entropy current have positive divergence: > 

0 . 

First, what is the entropy current? In ideal hydrody¬ 
namics 

= spu^ (29) 

where the entropy per baryon s depends on the equation 
of state P = P(p, u). Here we assume 

P = (7-I)u (30) 

and 

P = pQ. (31) 

The first law then implies ds = (du/u — 'ydpjp)l(^ ~ !)• 
One can show that, if r = proper time, 


Eckart frame eliminates another term related to mass 
diffusion). The factors of 1/0 are chosen for convenience. 
The ordering in the above expansion is ~ ^ 

e I. The first term spu^ is the leading order term 
O(e^) and is present even in the ideal case. The term 
oc is first order in a dissipative field 0(e), i.e. here 
the heat flux q>^. The terms oc q^'qa, II“^no ,/3 and 
are second order O(e^). Note that there is no viscosity 
contribution at first order in the above expansion, as is 
explained in section (I5.2|l . 

We now set ci = 0 to simplify the model further. The 
neglected term couples q and AP. The value of ci cannot 
be determined at this thermodynamic leve l and one has 
to res ort to kinetic theory (see section C in iBouras et al.l 
[2013). However, its choice does not affect the amount of 
entropy production. We remark further on the effect of 
Cl 7 ^ 0 at the end of the derivation. Thus we have 

= spu^ + 

where we have used H^'^H^j, = |AP^. Now evaluate 
First, 

yt.{spu^) = (36) 

where Tq+v is the sum of the conduction and viscosity 
terms in the stress-energy tensor. Then the conduction 
terms give 

UyV ^{u^q'" + q^^u'") =^q^. (37) 

where is the four-acceleration, and we 

have used the constraint u^q^ = 0. The viscosity terms 

give 

u^^,{-AP{m''-hn) = AP . 

In deriving the above, we have used the constraint 
= 0 uJfVvb^ = —b^b^'VvU^. Next, the first 
order term in q is 

(^) (39) 

the second order term in q is 


O = = (32) 

where TjOeai i® ideal gas stress-energy tensor. Com¬ 
bining this result with the continuity equation gives 

= 0. (33) 

In nonideal hydrodynamics one thinks of q^ and 
as small corrections to the ideal model. Expanding to 
second order in these small corrections, the most gen¬ 
eral possible entropy current subject to the constraints 
qf^ui^ = 0, H'/u'^ = 0 and H)) = 0 is: 

( 34 ) 

This is precisely what is done in Israel-Stewart theory, 
except that our bulk viscosity is 0 (and working in the 



biqdq q^ (blU^^\ 

e dr 2 0 ; 


(40) 


where dqjdr = and the second order term in AP 

is 


(v30 ) m dr 3 V 0 / 

(41) 

Assembling all the first-order terms in q from V^s^, 

-b . (42) 

The term proportional to has indeterminate sign, 

so we choose oi = 1. Gathering all the terms (first order 
+ second order) in q (and writing down only the q terms 
of {spu^)), we find 
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spu^ 




0 


20 


V ^0 

'W ~ ■© 



( biu‘^\ 

bibfj. 

dq 

1 0 ) 

0 

dr 


(43) 


Evidently this can be positive definite if the quantity in 
square brackets = PiQfj, and /3i > 0. Remarkably, by 
applying this condition we find an evolution equation for 

q- 


dq 

dr 


0 

h 


Piq- 


&'^(V ^0 + 0 a^) 
02 





(44) 

All that remains is to fix the constants bi and /3i. For 
/3i, we require that q asymptote to its non-relativistic 
value, —px^Qj where x is the conductive diffusivity 
and has dimensions of a length times a velocity. Then 
Pi = (PX©^)~^- For bi, we require that it be propor¬ 
tional to a relaxation timescale tr. Then bi = tr/ (px0). 
Gathering all together, our final evolution equation for 
the heat flux is 


dq q-qo qd f tr \ 


(45) 


J 


where 

go = -PX^'"(V;,0 -f 0a^) (46) 

is the classical Eckart heat flux projected onto the mag¬ 
netic field. The heat flux q relaxes to the first order 
(Eckart) heat flux go over a timescale tr. The addi¬ 
tional term on the right is a second order correction and 
is formally necessary to ensure the positivity of entropy 
production. The importance of this term in a full cal¬ 
culation can only be gauged by performing a calculation 
with and without this term. Notice that equation (1451) 
can be rewritten in the remarkably simple, scaled form 

dQ _ Q - Qo 

dr TR 

with Q = g(rfl/(xP 2 ))i/ 2 ^ 

Next, assemble all terms in depending on AP 

(and writing down only the AP terms of (spu^)) to 
find 


- ^AP'^uf^ 1 = AP 

OKJ 





fb2Ut^\ 

262 

dAP' 

/ 3 

K 0 J 

30 

dr 


(48) 


Evidently this will be positive definite if the quantity in 


square brackets = P 2 AP and /32 > 0. This provides an 
evolution equation for AP: 


dAP 

dr 




-/32AP-b-(6^S"V 





0AP^ (b2U^^\ 

2b2 ^ J ' 


(49) 


Again set the coefficients by requiring that 62 be propor¬ 
tional to a relaxation time tr and that AP asymptote to 
its classical non-relativistic limit = ?>pv[bb : Vv—^V-v) 
(see the Appendix), with v = kinematic viscosity, to find 
62 = tr/[2pv) and P 2 = {3pvQ)~^. Gathering all to¬ 
gether, the final evolution equation for AP is 


dAP _ AP- APo APd { '^r'\ 
dr Tr 2 dr ® \ piyP) ’ 


(50) 


where 


APo = 


(51) 


is a covariant generalization of the iBrasrinsl^ ()1965r i 


( psu>^ + -b>^ - ^^q^R 


which is positive definite, has the correct units, and 


model, in which collisions balance the forcing of 
anisotropy by the velocity field (see appendix A). As we 
already saw in (1451) where the heat flux relaxes to Eckart 
theory over a timescale tr, the pressure anisotropy AP 
also, relaxes to its first order value APq, with the second 
term on the right hand side of (l50l) being a higher order 
correction required for positivity of the entropy produc¬ 
tion. Notice that equation (1501) can also be written in 
simplffied, scaled form: 


dD D — Dq 

dT Tr 


(52) 


where D = AP{tr/ {pyP))^/"^ . 

It is also useful to gather the full, relativistic dissipa¬ 
tion function: 




Tr 

bpvQ 


AP'^u>^ 


g2 1AP2 
Xp0^ 3 upQ 


(53) 


reduces to the correct dissipation function in the non- 
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relativistic limit as we shall show in the next section. 
Note that the right hand side of (|53l) is a function of the 
full heat flux q and the pressure anisotropy AP, that are 
solved for using (l45ll and (l50l) respectively, and not the 
relaxed forms go and APq- The above value of entropy 
production is invariant under the choice of the cross cou¬ 
pling coefficient ci in (1551) . However, had we not set 
Cl = 0, the evolution equations for q (1^ and AP (l50|) 
would both have additional terms which couple q and 
AP to each other. 0 

The model given by equations (|45ll and (|5n|) is derived 
using precisely the same procedure as the Israel-Stewart 
model, but the complexity is greatly reduced because - 
thanks to the magnetic field - there are only two non¬ 
ideal degrees of freedom. We will also see below that 
the model is not subject to the linear instabilities of the 
isotropi c first order theory (where 6i = 62 = 0) discov¬ 
ered bv iHiscock fc LindblomI l)1985h . provided the damp- 
ing timescale tr is chosen appropriately, as found by 
iHiscock &: LindblomI (|1983ll and as we shall see in Sec. [3 

5. CONNECTION TO NON-RELATIVISTIC 

DISSIPATIVE THEORY 

In this section we compare the equations governing 
the entropy scalar s and the entropy current dSSl) 
to their non-relativistic counterparts. It is important 
to stress that the entropy scalar s and the entropy cur¬ 
rent are distinct quantities which obey separate evo¬ 
lution equations. Both of these have an analog in non- 
relativistic dissipative hydrodynamics. In equilibrium, 
the two quantities are related by s = —s^u^/p, but this 
is not true in general. This can be seen explicitly in (I35|) 
where there are second order differences (~ AP^) be¬ 
tween s and —s^u^/p. 

5.1. Entropy scalar s 

In non-relativistic dissipative hydrodynamics, 
the evolution equation for the entropy scalar s 
(|Landau k. Lifshitall987ll is 

Ds 

/90— — = —V • q — n : Vu (54) 


0 => The above equation is in¬ 

dependent of the model for q^ and H^"^, and only uses 
the form of the dissipative component of the stress- 

tensor. Thus, the equation is valid for both the first order 
Eckart theory as well as the second order Israel-Stewart 
theory. The difference between the non-relativistic equa¬ 
tion (|54[) and the relativistic equation (1571) . apart from 
the 3-derivatives transforming into covariant derivatives, 
is the presence of the term, which has no equivalent 
in the non-relativistic case. 

5 . 2 . Entropy current 

The equation for the divergence of the entropy current 
(1551) has been derived using a second order ansatz for 
the entropy current (1551) . However, the value of the en¬ 
tropy production is the same as in first order theories, 
two of which are the relativistic theory by Eckart, and 
the classical non-relativistic theory of dissipation. Below 
we show, starting from the evolution equation for the 
entropy scalar s from the previous subsection, that the 
entropy production in (1531) is true for both the relativistic 
and non-relativistic first order theories. 

• For the relativistic first order Eckart theory, we 
start with (l57)) and use the hrst order isotropic heat 
flux q^ = q(( = —pxh^^’^ (57^0 + 0a^), the first 
order shear tensor H^i, = Hq^;, = —2pi^{V^u'') = 
-pvh’^hP {WgUp -f VpUg - ^Ur^) 

(|Andersson fc Comer! 120071) . along with the 
relativistic continuity equation (SD to get 

QoQo/j. . IIo no/ii/ , . 

e)~pxe^ ^ 2 pvQ 

In the anisotropic case, we again start with (157)) - 
but now use the heat flux q'^ = qo and the shear 
stress H'"'^ = -APq where qo and 

APq are given by dJS]) and ([CT|) respectively, along 
with Q to get 


psu^ 


psw 


I' 


h^\ = 


Qo 


PA0^ 


APq^ 

SpiyQ 


(59) 


where D/Dt is the convective derivative d/dt + u-'V. We 
now derive the corresponding equation in relativistic dis¬ 
sipative hydrodynamics starting from (I32|) . which itself 
has been derived from the first law of thermodynamics. 
Proceeding to do so 

fj e 

p0_ = (55) 

= M,V^(q^u" + q"u'^+n'^") (56) 

= -V^q'^ - q'^a^ - (57) 

where we have used dsa and the constraint = 

^ The additional terms when the cross coupling coefficient 
Cl ^ 0 are of the form ^ = ... + (...)6'‘V,.AP + (...jAPV^S'' 
and = ... + (...)f)^‘Viiq-|-(...jqVfjS''. Therefore, q and AP are 
driven not only by gradients of the background thermodynamic 
quantities and b^b^'S /but by gradients of each other. 

Such terms become important in the collisionless limit. The cross 
coupling coefficients could potentially be used to derive a more 
accurate model where the dissipative fields are gn, gj^, APy and 
APj_, with the correct coupling between the fields. We leave this 
to future work. 


We see that the first order Eckart theory satishes 
(155)) with = psu^ + qo/0, and has the same 
amount of entropy production as the second order 
Israel-Stewart theory. However, in this case the en¬ 
tropy current is sourced by qo (I46|l and APq m 
which are directly related to gradients of thermo¬ 
dynamic quantities, as opposed to q and AP on 
the right hand side of (|55|l , which are solved for us¬ 
ing (l45)l and (150ll respectively. It can also be seen 
from the above why the viscous stress does 
not contribute at first order in the entropy expan¬ 
sion, since H^;^ = = —2pv{V^Uv) produces 

the right amount of dissipation /{2pvQ) 

on the right hand side, without any corresponding 
term in the entropy current on the left hand side 

of (1551) . 

• In the non-relativistic case, we start from (I54|l 
and use the classical heat flux q = —pxVO, 
the shear stress 11 = —2pi/(Vu}, which 

in component form is explicitly H^fe = 
-pi/ {dvi/dxk + dvk/dx^ - {2/3)6^kdvi/dxi) 
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(iLandau &: LifshitlllQSTD . and the non-relativistic 
continuity equation to get 


d{ps) 

dt 


+ V.(p.u + |) 


q- q 


n: n 

2puQ 


(60) 


The above is true for the anisotropic case as 
well, where we have the heat flux q = — pybb • 


V0, the shear stress 11 = ~APn 


6 b-1/3 


and 


APo = 66 ; Vt) — V • v/i (|Braginskiilll965l l. 

Clearly, (l60l) is the non-relativistic limit of (l58l) . 
where —>■ ps when u/c <C 1, since in this limit, 
vP 1 and ^ 0. 


Thus, the value of entropy production (right hand side of 
(l5^ ), which is second order in and , is the same for 
the entropy current expanded to first order (relativistic 
and non-relativistic), as well as to second order. How¬ 
ever, expanding the entropy current to third order leads 
to addit ional higher o rder terms on the right hand side 
of ([53]) (IE1 et a1.ll2ni?ill . 

6 . MODEL PARAMETERS: VISCOSITY, 
CONDUCTIVITY, AND RELAXATION TIME 

The parameters of our theory are t/j, z/, and y. The 
entropy production (|5^ can be interpreted in two ways, 
depending on whether the plasma is collisional/weakly 
collisional or collisionless. If the plasma is colli- 
sional/weakly collisional, it is the microscopic entropy 
production due to Coulomb scatterings and r/j is the 
Coulomb scattering time scale, i.e., the mean free time 
between particle-particle collisions. The transport coeffi¬ 
cients V, and X s-re set by this time scale tr- In particular, 
both V and x of order cItr in relativistic collisional 
kinetic theory: 

X = and v = il^c^hTR 

(61) 

where and ijj are constant dimensionless parameters 
and h = \ + ^u/[pc^) is the relativistic enthalpy. 

For a monoatomic ideal gas in the Chapman-Enskog 
theory </> = (15/4)')/;. In the non-relativistic Braginskii 
theo ry, where Coul omb interactions dominate, ~ A.lip 
(see iKulsrudI l2005l l. For a relativistic hard sphere gas 
(see the clear discussion of iCercignani fc Kremerll200If) 
z/ oc T for 0 = kT/{mc^) ^ I (m = molecular weight). 

In a collisionless plasma, because of the absence of 
Coulomb scatterings, the scattering time scale diverges 
tr ^ oo and the entropy does not increase. In this case, 
(ESI) is the increase in a coarse grained entropy and tr 
is the mean free time between wave-particle scatterings. 
This is because coarse graining the Vlasov equation leads 
to a Fokker-Planck equation with an effective collision 
operator, as is done in quasi-linear theory. Therefore, 
we still use the closures EH) but with a different inter¬ 
pretation of Tr compared to the collisional case. Wave- 
particle scattering is determined primarily by fluctua¬ 
tions in the electromagnetic field that have frequencies 
of order the cyclotro n frequency of t he particles of inter¬ 
est (ions, electrons) (|Kulsrudl l2005f) . Such fluctuations 
can either be produced by a cascade from larger scales 
or by velocity space instabilities that directly excite high 
frequency fluctuations. Unfortunately, the efficiency of 


wave-particle scattering by these processes is not fully 
understood. Moreover, the high frequency turbulent fluc¬ 
tuations that dominate scattering cannot be resolved in 
fluid (GRMHD) simulations, so subgrid models of the 
scattering rate are necessary. The formalism developed 
in the previous sections is sufficiently general that as the 
theoretical understanding of wave particle scattering de¬ 
velops, increasingly sophisticated models of tr can be 
implemented in our model. In particular, we stress that 
the relaxation time, viscosity, and conductivity in our 
model can be functions of the local plasma conditions 
(including, e.g., the local plasma /3, the amplitude of the 
turbulent fluctuations, etc.). Here we provide a rough 
guide to some of the key physics that motivates particu¬ 
lar choices for tr, v, and y. 

It is unclear to what extent magnetized turbulence 
in accretion disks (driven by the magnetorotational in¬ 
stability) produces significant wave-particle scattering. 
The energetically dominant component of the turbulence 
(associated with the slow and Alfven modes) does not 
produce efficient sc attering because it does not cascade 
to hig h frequencies l|Goldreich fc Sridharlll995l : IQuataertI 
1199811 . The fast mod e component can in prin ciple cascade 
to high frequencies (|Yan fc Lazarian| 120041) . but is very 
strongly damped in low-collisionality plasmas, which 
likely limits its importance. Absent signihcant scatter¬ 
ing by the high frequency tail of a turbulent cascade, the 
most important source of scattering is due to velocity 
space instabilities. These are instabilities in which tem¬ 
perature variations in different directions or local stream¬ 
ing velocities relative to the magnetic field frame pr ovide 
a sou rce of free energy to drive instabilities (e.g., iStixl 
[mi). Key examples include the firehose, mirror, and ion 
cyclotron instabilities for the ions and the electron fire¬ 
hose and whistler instabilities for the electrons (we show 
in (Elthat the firehose instability is captured by the fluid 
model developed in this paper; this is not true for reso¬ 
nant instabilities such as the ion cyclotron and whistler 
instabilities). The saturation of these instabilities is 
an active area of research with implications for galaxy 
cluster plasmas and the solar w i nd, as well as accre¬ 
tion disks fe.g.. iKunz et ^12014 IRiauelme et al.l[2015l : 
iSironi fc: Naravanl 120151 : iHellinger et al.l 1201511 . Note 
that for the ions in RIAFs, a non-relativistic the¬ 
ory {9i = kTi/{mi(?) < 1) is sufficient, since 6i ~ 
{H/r)'^{GM/{rc'^))^^^, and r > H = scale height. How¬ 
ever, for the electrons a relativistic theory is required 
since in the RIAF model > 1 and i ndeed relativistic 
electr ons are observed close to Sgr A* (jPoeleman et al.l 

[200l). 

Velocity space instabilities in low collisionality plasmas 
such as RIAFs are believed to arise because shearing, 
heating, expansion, or compression drives the plasma 
towards an unstable configuration, initiating the insta¬ 
bility. The scattering rate produced once the velocity 
space instabilities saturate is of the order of t he timescale 
on which the instabilities are driven (see IKunz et al.l 
12014 [Rlauelme et al.ll201^ . This balance between driv¬ 
ing and scattering timescales maintains the plasma near 
the marginal state for the instability. These consider¬ 
ations provide good motivation for one specihc choice 
of the relaxation time in low-collisionalty plasmas: tr 
of the order of the dynamical time r^, because the lat- 
















































Extended MHD 


9 


ter sets the characteristic timescale for heating, expan¬ 
sion, etc. For black hole accretion problems, this suggests 

TR^Td^ 

More specifically, however, velocity space instabilities 
will only set in if the free energy driving the insta¬ 
bility is sufficiently large. For example, the heat flux 
due to conduction cannot exceed the saturated value of 
q ~ where 3> is a dimensionless constant of or¬ 

der unity, but the electron whistler instability implies 
an even m ore stringent limit on the heat flux for high 
/3 plasmas (jPistinner fc Eichlerlll998l l. As a second ex¬ 
ample, a non-relativistic plasma is firehose (mirror) un¬ 
stable if AP/P < —2/j3 {AP/P > 1//3). Observa¬ 
tions of the solar wind show that the plasma pressure 
anisotropy obeys these constra i nts to reasonabl e accu¬ 
racy fe.g.. iHellinger et aril2006t iBale et aI1l2009ll . Thus 
a physical model of accretion disk viscosity should en¬ 
sure that |AP| does not exceed these bounds. Since, 
both viscosity and thermal conductivity are related to 
the relaxation time scale tr in (1611) . this implies a mod¬ 
ification of the thermal conductivity as well, when tr is 
decreased. 

One strategy for reducing tr in the presence of a small 
scale instability is to set 


Tr = Tdf f ——) (62) 

where x is some parameter (e.g. AP) that has a critical 
value Xcrit for instability. The function / is arbitrary 
but should have (1) /(O) = 1; (2) /'(O) = 0; (3) f{l) <C 
1. One function with the desired properties for / is the 
Fermi-Dirac distribution: 

= eO-i)A + l 

for y > 0, and f{y) = 1 for ?/ < 0. A is an ad¬ 
justable parameter that determines the width of the tran¬ 
sition to small Tr. As an example, saturated conduc¬ 
tion can be implemented by sending tr —?► Tdf{q/qcrit), 
where qcrit = is the maximum heat flux (and 

we have assumed a characteristic, unsaturated, relax¬ 
ation time of TdjEl. To consider several instabilities 
with instability threshold ratios yi, y 2 , ya---, one can set 
tr = Td/(yi)/(y2)/(y3)--- 

To understand the effect of (15^ and the associated 
reduction in the relaxation time near an instability 
threshold, consider an instability with threshold pressure 
anisotropy APcrit- As AP —>■ APcrit, the relaxation time 
is reduced. Because APq rx v rx tr, this reduces both 
APq and the time required for AP to relax to APq . Thus 
the plasma quickly falls below the instability threshold. 
Once below the threshold, tr becomes again, and AP 
relaxes to the large APq over a dynamical time scale. 
When the threshold is crossed again, (1621) takes effect 
and the cycle restarts. The net effect of this cycle is to 

I 


set up a feedback loop that results in q and AP being 
pinned to their values at an instability threshold, in a 
statistical sense. At every instant, the process generates 
entropy according to (15311 , which only involves q and AP 
and not qg and APq. 


7. LINEAR THEORY 

Here we address the following questions: (1) what are 
the characteristic speeds in our extended MHD model? 
These are needed to determine the Courant condition in 
explicit numerical evolution; ( 2 ) which choices of model 
parameters tr,x,v yield a stable model when the ini¬ 
tial state is an equilibrium? (3) under what conditions 
does one recover the firehose instability of an initially 
anisotropic state? (4) how large a heat flux can the model 
admit before it loses its stability and hyperbolicity? (5) 
what is the stability of the model in a frame not comov¬ 
ing with the fluid? and finally ( 6 ) under what conditions 
is the model causal? 

Consider an initially homogeneous magnetized fluid 
in Minkowski space. The initial state in the fluid 
frame has = { 1 , 0 , 0 , 0 }, p = po, u = ug, = 

{0,bsm9,0,bcos6), and q = AP = 0. The initial tem¬ 
perature is is Tg = Pg/pg = (7 — l)ug/pg. We per¬ 
turb around this initial state, e.g. q ^ 0 + Sq with 
Sq oc exp{ikx — iiot), linearize, and find the dispersion 
relation D{u!,k) = 0. 

It is worth first revisiting the linear theory for relativis¬ 
tic ideal MHD. As usual, the Alfven waves factor and the 
Alfven velocity is 


Va = 


pghg + 62 ■ 


(64) 


The slow and fast modes combine in a complicated, 
fourth-order dispersion relation. The special cases par¬ 
allel and perpendicular to the field give the sound speed 


.2 7^*0 

* Pohg 

and the fast magnetosonic speed 

2 _ 6 ^ + 7-Po 

Pohg + 6 ^ 


(65) 


( 66 ) 


respectively. 

A rigorous stability analysis would require a study of 
the general, ninth-order dispersion relation. Instead we 
consider the special cases of propagation parallel and per¬ 
pendicular to the magnetic field lines for conduction only 
and for viscosity only, and analyze the resulting disper¬ 
sion relations. 

ly = 0, k \\ B. The four Alfven modes are unaffected, 
as one might expect because they do not perturb the 
temperature. The remaining four modes arise from the 
entropy mode and the sound waves coupled to evolution 
of dq. The dispersion relation is: 


(w^ - clk'^)u}{ojTR + i) 



- lyky + Pk\uj^ - cik^) 


2,2n 7(7- l)c^ 


= 0 . 


(67) 


^ Saturated conduction is technically a collisionless effect and in 
this limit, the anisotropic pressure depends not only on the shear 


projected onto the field lines, but also on the gradients of the heat 
flux which we have ignored. See also footnote [2] 
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Evidently, the sound waves emerge in the ideal limit. in a numerical implementation of the model the Courant 

ly = 0, k J- B. We recover five zero frequency modes, condition must be adjusted accordingly, 

two magnetosonic modes, and one new mode: 

2 \-i 7.2. Stability 

■ (68) It is well kno wn that the Eckart first order conduc¬ 

tion model (see IMisner et al.l I1973L p. 567) is sub- 
X = 0, k \\ B. We recover four Alfvenic modes (prop- ject to a dramatic instability with growth rate oc l/y 

agating in each direction with two polarizations), a zero (|Hiscock fc Lindblomlll985fl . The Eckart theory is recov- 

frequency mode that asymptotes to the entropy mode, ered in the limit t/{ —>■ 0 in our model, so we expect that 

and modes that couple sound waves and AP: it will be unstable if we set th too small. 

1 / 4 \ 2^2 To make this idea more precise, consider the /c —>■ 0 

oj^+iu}'^ -w I -p- - - - -^ -h ) -i —— = 0. limit of dHT]). We find 

tr \ 3 tr{ 1 +xuo/poc‘‘) ) tr 

(69) / 

X = 0, k ± B. We recover a damped viscous mode u = —i Itr -1 1 (76) 

coupled to the two magnetosonic modes: ^ 7C / 



uj'^-\-iuj'^ ——oj 


vk'^PQC^ 


Tr V 3rfl (62 + pp c2 -I- yrtp ) 


7 2 2 

- k Vf^ 


jU2 2 

.ft, 


'm 


7.1. Characteristic Speeds 

The introduction of evolution equations for q and AP 
makes the equations hyperbolic and therefore also intro¬ 
duces new characteristic speeds. On dimensional grounds 
we expect a characteristic speed associated with conduc¬ 
tion oc x/tr and a characteristic speed associated 
with viscosity v\p oc vjTR. What are the constants of 
proportionality? 

First consider the pure viscosity case in the limit k —>■ 
oo. For parallel propagation 


2 7 2 
w = k 


C« + 


3 hi^TR 

and for perpendicular propagation 
1 V 


= k^ 




3 tr {1 + 'juo/po + 62/po) 
This motivates the definition 


«AP = o 


3 hoTR 


(71) 


(72) 


(73) 


as a characteristic viscous speed. In our closure model 
vap = Cg and in a numerical implementation the 

Courant condition needs to be adjusted accordingly. 

Next consider the pure conduction case in the limit 
k ^ oo. For parallel propagation the general solution 
for for the coupled entropy-conduction-sound wave is 
complicated and not much simpler than the dispersion 
relation itself, but the expression 


(cl + ± {c^, + (74) 

provides a reasonable first approximation (although it 
does not reveal that oc l/( 7(7 — l)c'^ — c^v^), which 
leads to slightly higher signal speeds when Cs,Vq ^ c). 
Here 

p2^(7-l)^ (75) 

^ tr 

is a characteristic conduction speed. For perpendicular 
propagation, the excitation of q is damped and nonprop¬ 
agating. In our closure model Vq = Cg-^/^(^T^^T) and so 


P where we have temporarily restored factors of c and = 

■ xPq/{po + luo/c^)- Evidently the model is unstable if 

tr < XclKlc^) (77) 

which is the same insta bility condition found by 
iHiscock fc LindblomI (|1985f) for an isotropic viscosity. 
For propagation perpendicular to the field, (I55)) . we re¬ 
cover the same stability condition. Setting y = (ftcIrR, 
instability requires (j) > for parallel or perpendic¬ 

ular propagation. 

In the comoving frame, the viscous modes are stable 
if AP = 0 in the initial conditions. All parallel and 
perpendicular modes are either neutral or damped. The 
stability for viscous modes in a non-comoving frame will 
be discussed in 17.51 

7.3. Firehose Instability 

In non-relativistic Braginskii theory, the inclusion of 
anisotropic conduction and viscosity opens up avenues 
for qualitatively new types of linear instabilities. For 
example, in the presence of anisotropic conduction, 
there are buoyancy instabilities driven b y temperature 
gradients rathe r than entropy gradients (jBalbusI 120001 : 
iQuataertI 1200811 . In the presence of anisotropic viscos¬ 
ity, the character of the magnetorotational instability 
can change, with viscous transport of angular momen¬ 
tum drivin g an instability even when magnetic tension is 
negligible ijOnataert et al.l [20021 : iBalbu^120041 1. And, fi¬ 
nally, a background pressure anisotropy can itself be sub¬ 
ject to a host of instabilities, at least one of which (the 
firehose instability) is well ca ptured by the EMHD-typ e 
models described here (e.g.. iSchekochihin et al.l [200911 . 
We expect that much of this rich physics will carry over 
to the relativistic theory described in this paper since it 
reduces to Braginskii theory in the non-relativistic limit. 
As a simple concrete example of this, we derive here the 
relativistic instability criterion for the firehose instability. 

Pressure anisotropy modihes the propagation speed of 
Alfven waves and, if AP is large enough, transforms 
them into unstable, non-propagating modes. Firehose 
instability has not appeared until now because our equi¬ 
librium has AP = 0. But we can recover the firehose 
instability by setting AP ^ 0 in the background state 
and taking r/j —>■ oo so that the pressure anisotropy can¬ 
not decay and the initial state is an equilibrium. 
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The resulting dispersion relation for parallel propaga¬ 
tion is: 


b'^ + Pa+ 7“o + \^P 

For AP = 0 we recover the Alfven waves. For AP > 0 
the propagation speed increases and becomes superlumi¬ 
nal if AP > |(po + 7 Mo). However, this would require 
AP/P > 37 /( 2(7 — 1)). For AP < 0 we recover the 
relativistic firehose instability criterion: 

+ AP < 0 ^ INSTABILITY. (79) 

This criteri on is c o nsiste nt with the relativistic kinetic 
criterion of iLerch^ (|1966f) . 

7.4. Instability for Large q/u 

iHiscock fc LindblomI (119881) (hereafter HL 88 ) have 
shown that the Israel-Stewart theory loses stability and 
hyperbolicity above a critical value of g/u ~ 0.08898 in 
the ultrarelativistic {u ^ p) limit. What does this imply 
for the EMHD model? 

Consider an initial state with a background heat flux 
q = go- We will perturb around this state to test stability. 
But to do this, we are obliged to take the limit —>■ oo so 
that the initial state with q = go becomes an equilibrium. 
It matters how this limit is taken: if we assume % oc r^, 
terms proportional to both t/j and x must be retained. 
Linearization around this equilibrium for modes with k || 
B yields a fourth-order dispersion relation for the coupled 
sound/entropy/conduction excitations. 

The stability of the q = go equilibrium depends on 7 , 
(j), q, p, and u in a complicated way. To compare with 
HL 88 we consider the ultrarelativistic limit with 7 = 4/3 
and u p; HL 88 also make a choice for the coefficient 
bi that is equivalent to setting cj) = 12/5. This simplifies 
the analysis, and one can show that the discriminant of 


the dispersion relation (which in this limit has real co¬ 
efficients) changes sign at qo/u = 0.08898, as in HL 88 . 
Indeed, one can show that HL 88 ’s A = 1 Israel-Stewart 
model is identical to ours if we restrict attention to mo¬ 
tion along the field. 

We can generalize HL 88 ’s analysis by allowing cj) to be 
different from 12/5. The critical qo/u = 2^/^/3 as ^ 0. 

Between ^ = 0 and (/ = 3 the critical qo/u is slightly 
above qo/u = (3 — (f>)/10, still in the ultrarelativistic 
limit. 

Is the instability a consequence of the second-order 
terms in the evolution equation for ql The higher or¬ 
der terms enter the linear theory only when go 7 ^ Oj so 
they have played no role until now. Turning these terms 
off and repeating the stability analysis, one finds the dis¬ 
criminant of the dispersion relation is positive definite. 
The instability is therefore seated in the second order 
terms. This is relevant for numerical implementations of 
the EMHD model, since it suggests that manipulating 
the higher order terms may improve stability, albeit at 
the cost of violating the second law. 

7.5. Stability in a non-comoving frame 

It is an interesting feature of the Israel-Stewart 
theory that in some cases perturbations comoving 
with the fluid appear stable, but even a small rel¬ 
ative velocity between the frame in which the per¬ 
turbation is define d and the fluid, makes t he per¬ 
turbation unstable (|Hiscock fc LindblomI B985D . Con¬ 
sider a background velocity u^ = (L, 0 , 0 ), with 

r = yr -|- the Lorentz factor, a background mag¬ 
netic field bfj, = {—bxU^/T,bx,0,0), and perturbations 
{Sp,Sua;,SAP) of the density, longitudinal velocity and 
pressure anisotropy. From the linearly perturbed equa¬ 
tions of mass conservation, evolution of the pressure 
anisotropy, and = 0 , we get (in the absence of 

heat conduction) 


Tdt{6p) = poU:cdt{Sux), (80) 

3poho{l + 2ul)dt{Su^) = 2u,r^dt{SAP)-3u,r^dt{Sp), (81) 

2poiyb2T^dtiSAP) = -T6AP + 2povu^dt5u:^. (82) 


Substituting for dtSux and dtSp, we get the evolution 
equation for SAP, 


dt{SAP) 


3{2po{ho - Bju/x -I- poho) SAP 

3b2{pohor'^ + po{ho - l)ul) - 2ul 2povT ' 

(83) 


In this case, the stability condition is 


b2 > 


_2u|_^ 

3{pohoT'^ + po[ho - l)ul)' 


(84) 


Note that 62 is the coefficient of the second order con¬ 
tribution to the entropy current (1551) due to shear stress 
and is equal to th/{ 2poi'). 


7.6. Causality 

For small r/j, i.e., for small bi and 62 , the model is no 
longer causal. To demonstrate this, we consider the prin¬ 


cipal part of the evolution equations at a point at which 
Ufx = (— 1 , 0 , 0 , 0 ), and 6 ^ = ( 0 , b^, 0 , 0 ), and compute the 
characteristic speeds of the system in the tetrad frame 
and along the direction . Writing the linearized sys¬ 
tem of equations as 


A*dt{SU) + A^d,{SU) + BSU (85) 


for perturbations SU, with matrices A^, B depending 
on the unperturbed variables Uo, then the characteristic 
speeds v along direction i are solutions of the equation 


det{v * A* — A*) = 0 . 


( 86 ) 
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In the absence of heat conduction, the relevant parts of 
the equations for (p, u, Ux, AP) are then 


dtP — pdxUx 

(87) 

dtu = -{u +P)dxUx + ... 

( 88 ) 

phdtUx = -dxP+^dxAP + ... 

(89) 

dtAP=^dxUx+ ... 

02 

(90) 


(the transverse components of the magnetic field and 
velocity still follow the standard dispersion relation for 
Alfven waves). The characteristic speeds of the system 
then depend on the equation of state P(p, u). For a sim¬ 
ple gamma-law fluid, P = (7 — 1 )m, we find the speeds 
(0, 0,±^/{2ph' -I- 3(7 — l)jub 2 )/{3phb2)). For small val¬ 
ues of & 2 , two of the speeds are thus greater than c = 1 . 
To avoid superluminal speeds, we need 


3(p-b 7 u (2 - 7 )) ■ 


(91) 


Rewriting in terms of the relaxation time scale r/j, the 
condition for causality is 


tr > 


4pR 

3(p-b7n(2-7))' 


(92) 


This is identical to what we would have inferred from 
equation (ITT]) by imposing uj/k < c. Regardless of 
our choice of frame, the theory is thus problematic for 
small 62 , i-e. small tr or large v/tr. These results 
are on ce more reminiscent of what iHiscock fc LindblomI 
(flMl found for an isotropic viscosity. In that case, 
they derived the stability condit ion &2 > 2 /( 3 / 9 /i). 
IHiscock fc LindblomI (119831 Il988bll also demonstrated 
that the stability of the theory implied that the system of 
evolution equations was hyperbolic and causal, and that 
if the system of equations was either not causal or not 
hyperbolic, it was unstable. 


8 . NONLINEAR THEORY 

In this section we give a brief but incomplete descrip¬ 
tion of the behavior of the model in shocks. 

Ideal hydrodynamics allows shocks and contact (en¬ 
tropy) discontinuities. The boundary conditions at the 
discontinuity are determined by continuity conditions on 
the fluxes of mass, momenta, and energy. 

Once viscosity is included in a non-relativistic fluid 
model the discontinuity is replaced by a sharp but con¬ 
tinuous transition of width Ax such that the diffu¬ 
sion timescale Ax^ jv is comparable to the transit time 
through the shock Axjv, so Ax ^ v/v. li v ^mfpCs 
and V ^ Cs then Ax ^ \mfp- Put differently, viscos¬ 
ity can propagate information upstream at character¬ 
istic speed ^ vjAx, which is supersonic if the struc¬ 
ture is narrow enough. Conduction produces an up- 
stream precursor but does not re move the discontinuity 
(|Mihalas fc Weibcl Mihaia^ll984[l . 

In Maxwell-Cattaneo type theories there is a new char¬ 
acteristic speed, oc (i^/ta)^/^; if 1 / ^ then this speed 
is comparable to Cg. If the upstream velocity exceeds 
this, the shock is still disc ontinuous. In Israel-St ewart 
theory, this was studied bv lOlson fc Hisco^ (II990D . See 


iBouras et al.l (120101) for a numerical study of shocks in 
this theory and a comparison to shock solutions obtained 
using kinetic theory. 

Linear analysis showed ()J7]) that our EMHD model 
contains a speed ~ ~ associated with 

thermal conduction and another speed ^ associ¬ 

ated with viscosity. One would then expect the model 
to have a discontinuity over a sufficiently strong shock. 
How, then, do AP and q change across the discontinuity? 

The variation of AP and q depend on the shock sub¬ 
structure, essentially because there are no continuity con¬ 
ditions, such as those present in ideal MHD, that can be 
used to navigate across the shock. To clarify this point, 
consider the conduction model governed by equation (I47|) 
in the frame of a strong (so that v upstream is greater 
than Vq) steady shock at x = 0, with > 0. The domi¬ 
nant terms near the shock are 



because both dxT and dxUx are large inside the shock. 
This equation has the form 

dxQ = F{p,Q,u^)5{x) (94) 

where 5 arises from dx acting on the discontinuities. If 
F were continuous through the shock we could integrate 
this equation across the jump and find a definite solu¬ 
tion for Q at X > 0. But F is changing discontinuously 
through the shock and hence the postshock Q depends 
on the substructure of the shock. 

EMHD shocks contain substructure on three scales. 

(1) On small scales the width of the shock is set by 
bulk viscosity. 

Because bulk viscosity is not included in the model, 
this smallest scale is unresolved and the evolution of 
q and AP across the shock is undetermined (but the 
change in, e.g., Q across the shock = AQ is limited, 
since if F in equation (IMl) is monotonic then AQ must 
lie between F in the downstream state and F in the up¬ 
stream state). In a numerical evolution it will be set by 
the numerical closure. 

(2) On intermediate scales the pressure anisotropy and 
heat flux relaxation produce a tail of width ^ u^tr. 

The separation of scales (1) and (2) is an artifact of 
the model since if tr is small inside the shock, (I) and 
( 2 ) can be comparable. 

(3) Far from the shock AP and g —>■ 0 and the shock 
obeys the ideal MHD jump conditions. 

The latter point implies that unless one is interested 
in shock substructure (and it would seem that EMHD is 
not the ideal model for collisionless shocks), the EMHD 
model may provide an adequate description of flow more 
than a few mean free paths from the shock. 

9. DISCUSSION AND CONCLUSION 

In this paper, we have derived an extended MHD 
model for a relativistic weakly collisional plasma incor¬ 
porating the effects of anisotropic conduction and vis¬ 
cosity. We are motivated by applications to low accre¬ 
tion rate black holes, but the model can also be applied 
in other contexts such as neutron star magnetospheres. 
This is the first in a series of publications. While this 
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paper deals with the derivation, linear theory and stabil¬ 
ity analysis of our model, a later publication will discuss 
a flexible new code which we have written to solve this 
model numerically, as well as possible future extensions 
of this model. As tests for the code, and to get insights 
into the model, we derive various analytic/semi-analytic 
solutions to this model. The code, as we ll as th e various 
tests will be described in iChandra et al.l (|2015ll . Finally, 
we shall apply this model and the code to study the dy¬ 
namics and structure of an accretion disk around a slowly 
accreting Kerr black hole. 

9.1. Summary of the Governing Equations 

The complete model is given by the usual continuity 
equation energy-momentum equations © and in¬ 
duction equation (TTKll . supplemented by a heat conduc¬ 
tion model (1251) and shear viscosity model (l28ll . together 
with evolution equations for the magnitude of the heat 
flux q (HSl) and the pressure anisotropy AP (1501) . 

9.2. Summary of the Formalism 

We deduced the form of the heat flux (|25|) and the 
viscous stress ()28[) by examining the symmetries of the 
system in the presence of a magnetic field and by appeal¬ 
ing to the small values of the gyroradius and the gyrope- 
riod in astrophysical plasmas. The evolution equations 
dUl) and (l50)) are then derived using the Israel-Stewart 
theory of dissipative hydrodynamics. We have expanded 
the entropy current up to second order in deviations from 
equilibrium and then enforced its divergence to be pos¬ 
itive. From this simple principle, we get equations for 
the heat flux and the pressure anisotropy in a magne¬ 
tized plasma. The derivation naturally shows that, as 
in non-relativistic plasmas, 1) the heat flux is driven by 
temperature gradients projected along the field lines and 
that 2) the pressure anisotropy is driven by a background 
shear projected along the field lines. Since our model is 
based on the Israel-Stewart theory, it conditionally satis- 
hes all the requirements of dissipative theories in general 
relativity: hyperbolicity, causality and stability. For ex¬ 
ample, in the ultrarelativistic limit {u ^ p), the model 
is hyperbolic as long as q/u < 0.089. 

In the limit where the relax ation time s cale t r —>■ 0, 
our model reduces to that of iBraeinsl^ (|I965ll . The 
hyperbolic equations for parallel heat flux and pressure 
anisotropy relax to forms ()46|) and 62) respectively, 
which are covariant generalizations of the Braginskii clo¬ 
sure. The pressure anisotropy in Braginskii theory arises 
due to a balance between the conservation of adiabatic 
invariants and pitch-angle scattering. In the appendix, 
we show using relativistic kinetic theory that this is true 
of our model as well. In the non-relativistic limit, the 
inclusion of anisotropic conduction leads to many new 
types of instabilities , such as the magnetothermal insta¬ 
bility (lBalbusll2000H and t he heat flux driven buoyancy 
instability (|Quataertll2008ll . both of which a re also mod i- 
fied by the inclusion of anisotropic viscosity (lKunzll20Tlll . 
Our relativistic model should also display all the instabil¬ 
ities that arise due to inclusion of anisotropic dissipative 
effects. As an example, we showed that our model repro¬ 
duces the correct threshold for the firehose instability. 

We relate the transport coefficients in our model (y, z/) 
to the effective scattering time scale t/j. A collision¬ 
less plasma is subject to a number of kinetic instabilities 


(mirror, firehose, ion cyclotron, electron whistler) that 
effectively regulate the pressure anisotropy and the heat 
flux. Our closure prescription incorporates the isotropiz- 
ing effects of these instabilities through the modification 
of Tfi, thus providing a convenient way to incorporate the 
various kinetic effects. The exact way in which tr should 
be modified will ultimately be answered by first princi- 
ple particle-in-cell calculations, which ar e well under way 
(iKunz et al.ll20l4 iRiauelme et al.l[20I5[) . 

9.3. Two-Fluid Effects and Observational Signatures 

The model is a one-fluid model of a plasma consisting 
of multiple species, here electrons and ions. The one-fluid 
matter stress tensor T|()atter (122 has been obtained by 
summing up stress tensors of each individual species in 
®. Therefore, the rest mass density p, the internal en¬ 
ergy u, the pressure P, the heat flux q^^ and the shear 
stress 11'^'^ in (1211) are all a sum of the respective quan¬ 
tities of both electrons and ions (using fv uf = u^), 
i.e. p = Pe + Pi, u = Ue -\- Ui, P = Pe Pi, q^ = q^ qf 
and = Wff -\- Tff’'. The temperature 0 in the model 
is a combination of both electron ©e and ion temper¬ 
atures 0i, i.e., 0 = {miQi me0e)/(mi -I- me) (since 
P = Pe Pi), where m-e and mi are the electron and ion 
masses, 0e = kTe/{me<?), and 0^ = kTi/{miC^). Intro¬ 
duction of a single heat flux q^^ and a single shear stress 
for a system with multiple species is inconsistent, 
unless we make further approximations. 

The ratio of the relaxed electron to ion heat flux is 
qeo/qio ^ {PelPi)'^{Teln){mJme){LilLe), where Le and 
Li are the characteristic scales for the variation of the 
electron and ion temperatures and Te and Ti are the 
electron and ion relaxation time scales respectively. As¬ 
suming Le ~ Li and Te Ti, our one-fluid model is 
a reasonable description of the overall dynamics when 
Pe/Pi \fmfjmi. In this limit, the contribution to 
the total heat flux is dominated by the ion heat flux 
q^ ~ q^, and the single temperature 0 in our one-fluid 
model is the ion temperature 0 ~ 0^. Similarly, the 
ratio of the relaxed electron to ion pressure anisotropy 
is APeo/APio ^ (PelPi){reIn) < I, when Pe/Pi < I 
(again assuming Te ^ n). In this limit, the total pres¬ 
sure anisotropy is dominated by the ions and thus we 
have for the total shear stress 11^'^ ^ n(“^. Note that 
in the collisional (Braginskii) regime, Pe — Pi, and 
Te/n ^ \fmfjmi, where Te and Ti are the Coulomb scat¬ 
tering time scales for electrons and ions. The dominant 
heat flux is now due to electrons geo /<?io ~ \Jrni/me ^ 1, 
while the dominant pressure anisotropy is still due to ions 
APeo/APiQ ~ yjmfjmi ^ 1. Since in this collisional 
limit, Te — Ti, evolving a single temperature 0 20^, 

as we do in this paper is appropriate. 

Generalizing this one-fluid model to a two-fluid model 
that accounts for both electron and ion dissipation can 
be done in a number of ways. One approach is to use 
separate conservation equations for both the electrons 
and ions, into which our model can be incorporated in 
a straightforward manner. This however includes kinetic 
length scales in the system, which cannot be resolved in 
a global accretion disk simulation. Another approach 
is to work within the MHD ordering, which only re- 
quires one additiona l variable, the electron temperature. 
iRessler et al.l l)2015ll have developed such a model by us- 
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ing a separate electron entropy equation into which they 
incorporate electron conduction using a reduced version 
of the anisotropic conduction equation that we have de¬ 
rived here. The emission from the plasma is dominated 
by the electrons and hence the electron thermodynami¬ 
cal quantities are crucial to predict observables such as 
spectra, images and light curves. 

The non-ideal effects modeled here may have a number 
of implications for the dynamics and observational prop¬ 
erties of slowly accreting black holes: 1) conduction can 
redistribute energy spatially, changing where the emis¬ 
sion comes from and potentially changing the dynamics 
by increasing the pressure at large latitude/radii and 2) 
viscosity can provide an additional source of transport 
and heating (over and above the Maxwell and Reynolds 
stress), potentially modifying M and the thermodynam¬ 
ics of the plasma - the latter being relevant for the emis¬ 
sion. 

9.4. Connection to Other Models 

The model we have presented is the simplest in a hier¬ 
archy of possible models incorporating non-ideal effects 
in magnetized plasmas. In particular, it only includes 
a single parallel heat flux and the parallel 

and perpendicular pressures are controlled by a single 
variable, the pressure anisotropy AP. The derivation 
was based on thermodynamic principles. On the other 
hand, models of collisionless plasmas are usually derived 
by taking velocity space moments of the Vlasov equa¬ 
tion. Earl y examples of re l ativistic models der i ved in this 
way are bviGedalid (11991 11 . iTsikarishvili et al.l(|1992ll and 
iTsikarishvili et ahl (|1994f) . They have independent vari¬ 
ables for the parallel (P||) and perpendicular (P_l) pres¬ 
sures. However, their stress tensors lack a heat flux and 
their models reduce in t he non-relativistic limit to the 
well-known C GL closure (lOiew et al.lll956ll . A more re¬ 
cent model bv lTenBarge et al.l (|2008ll has evolution equa¬ 
tions for P|| and P_l and these couple individually to sep¬ 
arate heat fluxes ^ 6^V^P|| and q± ^ 6^V^Px. In 
the non-relativistic regime, when the mean free path is 
small compared to the s ystem size, such a m odel recov¬ 
ers the Braginskii limit (|Snvder et al.l[T99^ . Thus, we 
expect our model to be formally applicable only in this 
regime. Also, since we have not derived our model using 
kinetic theory, we do not expect that it will reproduce 
all of the kinetic linear modes; for example, our model 
lacks the mirror instability. However, a nice feature of 
our model is that it includes collisions whereas the pre¬ 
vious relativistic models for collisionless plasmas do not. 
To include collisions in the moment formalism that the 
previous models used, one has to include a collision oper¬ 
ator in the Vlasov equation and take its moments. Our 


model naturally includes collisions because it is an ex¬ 
pansion up to second order around thermal equilibrium. 
The relaxation time scale tr can be interpreted as the 
relaxation time scale in the BGK collision operator. Its 
presence is convenient since it makes it easy to incorpo¬ 
rate subgrid models of the saturation of kinetic plasma 
instabilities as an effective collisionality. 

Our model captures the leading order effects of heat 
transport and pressure anisotropy while still being rel¬ 
atively simple. More sophisticated models have been 
derived from Isr ael-Stewart theory in the presence of 
a magnetic field (|Huang et al.ll2010ll . albeit in a differ¬ 
ent context, for use in strange quark stars with strong 
magnetic fields. The Israel-Stewart theory itself is a spe¬ 
cific instance of a class of theories derived from Extended 
Thermodynamics, which are based on an entropy prin- 
ciple, such as th e one we have used in this paper. See 
iJou et all (|1988ll for a review. While this formalism gives 
us the evolution equations for the dissipative fluxes, one 
needs to eventually resort to kinetic theory in order to 
compute the transport coefficient s. To perform this com¬ 
putation, [Israel &: Stewartl (1197911 have used an approach 
similar to that of IGradI (119491) . They expand the distri¬ 
bution function in momentum space polynomials around 
an equilibrium, i.e., / = /o(l -I- -I- where 

/o is the equilibrium distribution function. This ansatz 
for the distribution function is then used along with the 
second moment of Boltzmann equation to compute the 
transport coefficients, which come out in terms of mo¬ 
ments of /q. Note that the second moment of the Boltz¬ 
mann equation is one moment higher than the divergence 
of the stress-tensor and, it is at this level of the moment 
hierarchy where the collision operator makes a non-zero 
contribution. 
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APPENDIX 

ADIABATIC INVARIANTS AND PRESSURE ANISOTROPY 

Equ ation 1511 in the main text for APq is the natural general relativistic generalization of the well-known iBragins^ 
(j 19651) relation between pressure anisotropy, viscosity, and shear stress in a non-relativistic plasma. In the non- 
relativistic limit, this relationship also has a simple interpretation in terms of a balance between (1) the rate of 
generation of pressure anisotropy by adiabatic invariance of the magnetic moment in a slowly varying magnetic field 
and (2) the isotropization in velocity space by pitch-angle angle scattering at a rate Vg. We briefly review this non- 
relativistic result and then discuss its generalization to relativistic kinetic theory of a magnetized plasma. This provides 
a useful physical interpretation of the thermodynamic derivation of the equilibrium pressure anisotropy (eq. 1511) in 
































Extended MHD 


15 


the main text. We take k = m = c = 1 throughout this Appendix. 


Non-relativistic Kinetic Theory 

The viscous stress tensor in a magnetized no n-relativistic colli sional plasma in which the cyclotron frequency is much 
larger than the collision frequency is given by (jBraginskiilll965[ l 


Equation [AT] can also be written as 


n = — 3 pz/ 

bb : 'Vv - 

V • V 


3 



r— I] 


n 

= -AP 

“-3 



where the pressure anisotropy in the non-relativistic limit is given by 


AP = Spi' 


bb : 'Vv — 


V • V 
3 


(Al) 

(A2) 


(A3) 


AP in equations IA2I and IA3I is defined as AP = P± — Py, with directions defined by the local magnetic field in the 
plasma. Note that since v c'^/uj where tv is the pitch angle scattering rate and Cg is the sound speed, equation I A3I 
is equivalent to 


AT 


dt 


■In 


p3 


(A4) 


where we have used the i ndu ction equation and mass conservation to rewrite the right-hand-side of equation I A41 The 
left hand side of equation |A4] is the rate at which scattering at rate w decreases the pressure anisotropy. The right hand 
side of equation IA4I is the rate at which adiabatic invariance of T±/B (magnetic moment) and T^B "^/(the ‘bounce’ 
invariant) change the pressure anisotropy. The assumption of collisional theory is that these two effects approximately 
balance each other. 


Relativistic Kinetic Theory 

The thermodynamic derivation in 0 (i^i particular eq. El]) shows that results very similar to equations IA3I and IA4I 
relate pressure anisotropy and viscosity in GR. To understand the microscopic origin of these results, it is instructiv e 
to consider the relativistically correct first and second adiabatic invariants for an individual particle (|Sturrocklll994f) : 

tP" V\\b 

= constant - = constant (A5) 

b p 

where pj_ and py are the particle momenta (not pressure!) along and perpendicular to the magnetic field direction 
defined in the fluid frame. What does this imply for the relation between AP and the variation of b and pi 
Start by changing momentum space coordinates to the adiabatic invariants 


• _ pi 

j± = — 
m 

III 

(A 6 ) 

where r = p/p^ and m = b/b^. The initial {r = m = 

1 ) Maxwellian distribution function is 


dn 

^ -r /0 

(A7) 


40^2(1/0) 

where 



r = 

(1 + 

(A 8 ) 


and K 2 is the modified Bessel function of the second kind. This distribution is invariant under slow (adiabatic) changes 
in b and p. 

Next, directly evaluate components of the pressure tensor using the kinetic theory definition (HD in the limit that 
r and m are close to 1 : 


AP = Pi - P|| 





(A9) 


where the factor of ^ arises because p\ = p^+Py if 6 is aligned with z; here p* = T in the fluid frame. Notice that p_L 
and p|| depend on r and m but the distribution function does not. Then 

dAP 
dr 


dAP dm dAP dr 


ad 


dm dr 


dr dr 


(AlO) 
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is the change in AP due to adiabatic deformation of the distribution function. If the anisotropy is small the derivative 
can be evaluated at b = m = 1. 

The integrals needed for dAP/dm and dAP/dr can be evaluated analytically in the relativistic and non-relativistic 
limits. 

In the non-relativistic limit T —l + + and dn/dj±dj\\ oc 0“^/^exp(—(j_L-|-j||^))/(20). A few strengthening 

integrals later, one finds 


1 

P[l^ 


ad 


„ dm „ dr 


(All) 


where P = n0. Since the derivatives are evaluated at m = r = 1, dm/dr = dlnm/dr, etc., and using the definition 


of 1 


find 


1 /dAP\ 


P 




ad 


^ 1 


(A12) 


as expected from the discussion in the preceding subsection. 

In the ultrarelativistic limit T — {p^ = {j± + at 6 = m = 1, and dn/dj±dj^^ oc 0“^exp(—r/0). The 

integrals give 

1 / ^ \ TD\ 10 C 

(AI3) 
(AM) 


or 


dAP^ 

\ 12 dm 

dr 2 

Kd 5 dr 

dAP" 

) pti'" 

dr y 


which differs by a factor of 4/5 from the non-relativistic limit. 

The rate of decay of the pressure anisotropy due to scattering is 


1 fdAP\ 


P 


AP 




= —UJ 


which defines the scattering rate w, so the total change 

1 /'dAP\ 1 /dAP\ 


P \ dr 


') P\ dr ) 


ad 


p ’ 

1 fdAP\ 


In equilibrium scattering balances adiabatic forcing and, for the ultrarelativistic limit, the result is 


AP 4 d 


P 


= —In ^ 


5 dr 






(A15) 


(A16) 


(A17) 


We can do just a little bit more by evaluating the right hand side using the induction equation and the continuity 
equation, together with = 0. Along the way, 


and 


so that 


z dr 


— \np= 


^ In ( 5 1 = 3 ( 


dr 


Then in equilibrium 


M0, 


AP = 3p[--] ( --V, 

0 UJ 


This is consistent with the non-relativistic derivation in equation IA31 and, if we set 

40 2 

v = -oc rc„ 

5 u} 

with the relativistic thermodynamic derivation in 21 of fh® main text (in particular, equation [51] and surrounding 
results). 


(A18) 

(A19) 

(A20) 

(A21) 

(A22) 
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